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ABSTRACT 

We consider an analytic model of cosmic star formation which incorporates supernova 
feedback, gas accretion and enriched outflows, reproducing the history of cosmic star 
formation, metallicity, supernovae type II rates and the fraction of baryons allocated 
to structures. We present a new statistical treatment of the available observational 
data on the star formation rate and metallicity that accounts for the presence of pos- 
sible systematics. We then employ a Bayesian Markov Chain Monte Carlo method to 
compare the predictions of our model with observations and derive constraints on the 
7 free parameters of the model. We find that the dust correction scheme one chooses to 
adopt for the star formation data is critical in determining which scenario is favoured 
between a hierarchical star formation model, where star formation is prolonged by 
accretion, infall and merging, and a monolithic scenario, where star formation is rapid 
and efficient. We distinguish between these modes by defining a characteristic mini- 
mum mass, M > 10 11 M Q , in our fiducial model, for early type galaxies where star 
formation occurs efficiently Our results indicate that the hierarchical star formation 
model can achieve better agreement with the data, but that this requires a high effi- 
ciency of supernova-driven outflows. In a monolithic model, our analysis points to the 
need for a mechanism that drives metal-poor winds, perhaps in the form of supermas- 
sive black hole-induced outflows. Furthermore, the relative absence of star formation 
beyond z ~ 5 in the monolithic scenario requires an alternative mechanism to dwarf 
galaxies for reionizing the universe at z ~ 11, as required by observations of the mi- 
crowave background. While the monolithic scenario is less favoured in terms of its 
quality-of-fit, it cannot yet be excluded. 
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1 INTRODUCTION 

Massive galactic spheroids form either by hierarchical build- 
up or monolithically. The former scenario is favoured by the 
theory of, and evidence for, cold dark matter (CDM), the 
latter by some observations, most notably the evidence for 
down-sizing in stellar mass and on chemical evolution time- 
scales from the [a/Fe] enhancement at high spheroid masses 
(Worthey et al. (1992)). The conventional view of hierar- 
chical build-up fits in with semi-analytical galaxy formation 
simulations provided that suitable feedback models are pre- 
scribed (Croton et al (2006), Bower et al. (2006)). Gas 
cooling drives star formation below a critical halo mass and 
AGN quenching of star formation occurs at higher masses 
where the cooling is inefficient. This model reproduces the 
shape of the galaxy luminosity function at the high mass 
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end, and also produces old, red massive galaxies. Not yet ad- 
dressed in these studies is the issue of whether these massive 
galaxies can be formed sufficiently rapidly, the resolution of 
which may lend support to a monolithic formation scenario. 
Major mergers are a plausible ingredient of a monolithic 
model, for triggering rapid star formation. 

Not all studies agree on the role of major mergers. Not 
all luminous starbursts are triggered by such mergers. In- 
deed, in the case of disk galaxy formation, their role is 
likely to be small. There is evidence that major mergers 
are not exclusively responsible for the dominant star forma- 
tion episodes in starbursts forming stars at up to 200 Mo/yr 
where thick disks are seen (e.g. Hammer et al. (2005)). Such 
high star formation rates, also found in disks at z ~ 2 by 
Forster-Schreiber et al. (2006), may favour an interpreta- 
tion in terms of bulge formation since the disks are already 
present. The observed major merger rate is small out to 
z ~ 1.2 despite the strong increase in the comoving star for- 
mation rate (Lotz et al. (2006)). Of course, ULIRGs, with 
star formation rates of £ 5OOM0 /yr, are examples of ma- 
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jor merger-induced star formation. However, for the bulk of 
star formation in the early universe, there is little evidence 
to suggest that major mergers play a significant role. This 
motivates a hierarchical picture of prolonged minor mergers 
to build up mass. 

One can ask more generally whether minor mergers or 
gas cloud accretion are responsible for the gas supply. With 
regard to minor mergers, the answer seems to be negative, 
because the time-scale for supplying the gas would be many 
dynamical time-scales. Yet there are two persuasive argu- 
ments for a short star formation time-scale in massive galax- 
ies. SED modelling suggests down-sizing, the most massive 
galaxies forming first and with star formation time-scales 
of order a dynamical time. The [a/Fe] enhancement with 
increasing spheroid mass independently reinforces this con- 
clusion, since most of the star formation must have occurred 
before dilution of the star-forming gas by supernovae type la 
ejecta occurred. This result is not, or at least not yet, found 
in semi-analytical modelling of massive galaxy formation. 
It provides a strong argument for monolithic formation of 
massive galaxies, although one cannot of course exclude the 
possibility that hierarchical formation models will reproduce 
a similar result once more complex star formation rules are 
introduced. Indeed, simple spherically symmetric accretion 
models have been proposed (cf. Birnboim et al. (2007)) 
that allow the gas to accumulate in an isothermal halo at- 
mosphere prior to an accretion-triggered burst of star for- 
mation. In essence, this approach introduces a monolothic 
formation model in combination with hierarchical gas ac- 
cumulation. Unfortunately there is little evidence today for 
such gas-rich halos, suggesting that this process, if impor- 
tant, would only have been influential at very early epochs. 

To save hierarchical models, one needs an efficient way 
of converting gas into stars. The typical gas-to-star e- 
folding time in nearby well-studied sites of star formation, 
namely spiral galaxies, is several Gyr. By contrast, at a red- 
shift of 2 — 3, the star formation time-scale is as short as 0.2 
Gyr from both spectrophotometric SED filling (Maraston et 
al. (2006)) and [a/Fe] analyses (Thomas et al. (2005)). 
The AGN phenomenon is invoked for stopping star forma- 
tion by quenching the gas supply. This may keep massive 
galaxies red, but may not be enough to produce massive 
galaxies sufficiently early. One resolution may be to invoke 
positive feedback from AGN-driven outflows that overpres- 
sure protogalactic clouds and trigger star formation on a 
rapid time-scale (Silk (2005)). Whether or not this partic- 
ular solution turns out to resolve the dilemma (if indeed the 
problem persists in the perspective of improved star forma- 
tion modelling) is not the issue we focus on in this paper. 
Rather we revisit the case for hierarchical versus monolithic 
galaxy formation in terms of accounting for the data on star 
formation rate and chemical evolution. 

In this paper we apply a novel statistical approach 
that rigorously treats the key parameters in semi-analytical 
galaxy formation theory, with emphasis on reproducing the 
cosmic star formation and chemical evolution histories. Nu- 
merical modelling via semi-analytical techniques of a large 
box of the universe takes up so much computer time and 
memory that it is impossible to test the robustness of the 
results. Here we focus on probing the key parameter space by 
means of analytical techniques combined with appropriate 
binning of the full data sets. We find that the characteristic 



minimum mass for the building blocks of massive galaxies 
plays a central role. Specifically, the dust correction scheme 
one chooses to adopt for the cosmic star formation history 
data is one of the most critical factors in determining the 
balance of evidence in support of a hierarchical star forma- 
tion model as opposed to a monolithic scenario, where star 
formation happens predominantly in massive spheroids. Our 
results indicate that the hierarchical star formation model 
can achieve better agreement with the data, but that this 
requires a high efficiency of supernova-driven outflows. In a 
monolithic model, our analysis points to the need for a mech- 
anism that drives metal-poor winds, perhaps in the form of 
supermassive black hole-induced outflows. While the mono- 
lithic scenario is less favoured in terms of its quality-of-fit, 
it cannot yet be excluded. 

This paper is organized as follows: our star formation 
model is introduced in section [2] and the dependence of the 
observable quantities on the free parameters of the model 
explored in section [3] We then describe the data employed 
and our statistical procedure in section Our results are 
presented in section [3] and our conclusions discussed in sec- 
tion [6] Appendix [X] gives details of our binning procedure 
for star formation rate and metallicity data which accounts 
for undetected systematics. 



2 STAR FORMATION MODEL 

In this section we present a physical model of the cosmic 
star formation incorporating supernova feedback, gas ac- 
cretion and enriched outflows. Our model builds upon the 
model described in Daigne et al. (2005), and more specifi- 
cally their "Model 0" . The main difference at the model level 
is our choice of using a different chemical evolution model. 
In fact, we adopt the instantaneous recycling approximation 
whereas in Daigne et al. (2005) the metallicity is computed 
under the delayed enrichment approximation. There are sev- 
eral differences in terms of the statistical treatment of the 
data and the fitting procedure, that in this work are sig- 
nificantly more sophisticated, as explained in section U] and 
Appendix fX] 



2.1 Governing equations 

The description of baryons in the Universe and the pro- 
cesses that define the evolution of the baryonic mass are of 
fundamental importance for our model. Following Daigne et 
al. (2005) , we employ three baryon reservoirs in the model, 
encompassing the interstellar medium (gas), the mass in 
stars and the intergalactic medium (IGM). We denote by 
Mg as the mass of gas, by M* the mass in stars and by 
Mstruct = Mgas + M* the total mass in collapsed structures. 
The IGM and the structures exchange mass through accre- 
tion and outflow, while the interaction between stars and 
gas is governed by star formation and ejection of enriched 
gas. In the instantaneous recycling approximation adopted 
here, the accretion rate of the mass in stars is simply equal 
to the star formation rate (SFR), ^(t), i.e. 
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We then have the following set of differential equations gov- 
erning the evolution of the mass in the three reservoirs: 



dt 

It 



dM si 



dt 

a b (t) — o(t) 

dM s truct dM* 



dt 



dt 



dt 



(2) 
(3) 
(4) 



In the above equations, a,b(t) is the rate of baryon accre- 
tion while o(t) is the rate of baryon outflow. The latter 
includes two terms, accounting for winds powered by stel- 
lar explosions and supernova ejecta. We neglect supernova 
ejecta since this effect was found to be subdominant (Daigne 
et al. (2004)). The relation between physical time t and red- 
shift z is given by 



dt 
dz 



9.78ft _1 Gyr 



(l + ^) X /fiA+On(l+2) 3 



(•») 



where we have assumed a flat Universe. In the following, we 
will use Eq. ([5} with parameters fixed to the values of the 
ACDM concordance model, i.e. a matter density parameter 
fi m = 0.27 and a cosmological constant energy density 
Q.A = 0.73 (both in units of the critical energy density 
of the Universe), and we will take for the Hubble constant 
H = 100ft km/sec/Mpc = 71. 



2.2 Accretion 

We adopt the hierarchical scenario of structure formation, 
where small structures are formed first. At redshift z, the 
comoving density of dark matter halos in the mass range 
[M, M + dM] is f ps (M, z)dM, normalized in such a way that 



/oo 
dMMf pB (M,z) = Pdm, 
o 



(6) 



where pdm is the comoving dark matter density. The dis- 
tribution function of halos f pa (M,z) is computed using the 
method described in Jenkins et al. (2001) using code pro- 
vided by A. Jenkins. It follows the standard theory (Press & 
Schechter (1974)) including the modification of Sheth & 
Tormen (1999). We assume a primordial power spectrum of 
fluctuation with a power law index ns = 1 and the fitting 
formula to the exact transfer function for non-baryonic cold 
dark matter given by Bond &Efstathiou (1999). For the 
rms amplitude we adopt a value erg = 0.9 for mass density 
fluctuations in a sphere of radius 8ft _1 Mpc. 

Using the above expressions for the distribution func- 
tion of dark matter halos, we can calculate the fraction of 
baryons at redshift z that is allocated to structures, by as- 
suming that the baryon density is proportional to the dark 
matter density, with a proportionality factor given by the 
ratio of visible to dark matter density - in other words, we 
assume that light traces matter with no bias. The fraction 
of baryons in star-forming structures at redshift z is then 
given by 



/bar (2) 



f^. n dMMf pB (M,z) 
ffdMMf pa {M,z) 



(7) 



where M m i n is a free parameter controlling the minimum 
mass (in units of solar masses) of the collapsed structures 
where star formation can occur. 



The accretion rate is then given by (Daigne et al. 
(2005)): 



(Lb(t) — Qb 



3#o 
SttG 



I rf/ba 



dz 



(8) 



Given a value of M m i n (that we adopt as a free parame- 
ter, see below), we fix the redshift at which star formation 
begins, Zj n i t by the requirement that /bar(zinit) = 0.01. In 
other words, the first stars form in collapsed haloes of mass 
larger than M m i n when the fraction of baryons allocated to 
such structures is more than 1%. We adopt a fixed baryonic 
density parameter of fib = 0.044 (from the posterior mean 
of WMAP 3-years data combined with all other datasets, 
Spergel et al. (2006)). 



2.3 Outflow 

The adopted stellar initial mass function (IMF) is of the 
form 



$(m) = B 



Mr. 



-(1+x) 



for mi < m < m u (9) 



where the normalization constant B is fixed by the require- 
ment that 



m&(m)dm — M t 



(10) 



For the limits of integration we fix mi = O.IMq and 
m u = 100 M© (Pagel (1997)). Therefore the only param- 
eter needed to define the IMF is its power-law index, x. 
The quantity x is used as a free parameter in this model. 

We model the outflow powered by stellar explosions as 
follows: 



o(t) 
where 



2e /-IOOMq 
v 2sc( z ) Jm 



dm$(m)t>(t~T s (m))E kin (m) (11) 



mo — max(8M0, rrid(t)) 



(12) 



and <I>(m) is the IMF defined above, r s (m) is the lifetime of 
a star of mass m and md(t) is the mass of stars that die at 
age t. Furthermore, Eki n (m) is the kinetic energy released 
by the explosion of a star of mass m, that we take to be a 
fixed constant independent of mass, _Ekin(w) = 10 51 ergs (a 
mass-dependence could easily be taken into account). The 
free parameter e controls the fraction of the kinetic energy of 
supernovae that is available to power the winds, and Vg SC (z) 
is the mean square of the escape velocity of structures at 
redshift z. 

In order to compute the stellar lifetime t s (m) we assume 
it to be equal to the time that a star of mass m spends on 
the main sequence. So the age of a star of mass m is given 
by 



,(m) = (m/M Q ) 2 -H Q 



(13) 



where <q is the total time that a star of mass M = Mq 
will spend on the main sequence and we adopt a value 4q = 
9 Gyr. To compute m d (t) in Eq. (fT2")l we solve Eq. JT3J) for 
m, thereby obtaining the mass of stars md(t) that die at age 
t. 

The escape velocity is obtained by assuming virialized 
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halos and averaging over the distribution function, thus ob- 
taining: 

Jg . dMMf ps (M,z)(2GM/R(M)) 
Vcsc{Z) ~ JZ mm dMM,f ps (M,z) (14) 

where R(M) is the radius of a dark matter halo of mass M 
given by the following expression: 

/ 3^ \ 1 / 3 

R{M) = (,178^0(^(1 + ^+^)477 ) ' (15) 

where p c is the critical density of the universe to- 
day.The factor 178 is the overdensity (relative to the 
critical density) at virialization for an Einstein-de Sitter 
model (Coles & Lucchin (1995)). 

2.4 Star formation and supernova rate 

Following Daigne et al. (2005), we adopt an exponentially 
decreasing SFR: 

*(t) = t/Mrtruot(t) exp(-(i - t ini t)M, (16) 

where tjnit is the time corresponding to the redshift ziuit 
when star formation starts (as defined above), r is a char- 
acteristic time scale that we take as a free parameter and 
v is a normalization parameter (with dimensions of inverse 
time). 

The supernova rate (SNR) is strongly linked to the star 
formation rate because of the short lifetime of massive pro- 
genitors with M > 8Mq . We can therefore assume that core 
collapse supernovae (SNe) are strongly correlated with in- 
stantaneous SFR, and the supernovae type II rate $,„(t) is 
given by 

/•m u 

*sn(t)= / $(m)*(t-r,(m))dm. (17) 
JsM e 

2.5 Chemical evolution model 

Chemical evolution is included in the model using the in- 
stantaneous recycling approximation, i.e. we assume that all 
processes involving stellar evolution, nucleosynthesis and re- 
cycling take place instantaneously on the timescale of galac- 
tic evolution. The equation of galactic chemical evolution is 
(Pagel (1997)) 

Eg ^ = i 9 ® + {Zf ~ z (*))° b W - fa - WM*). ( 18 ) 

where E g is the density of the gas (in units of A/q/Mpc 3 ), 
r) is a multiple of the nucleosynthetic yield that parame- 
terises the metallicity Z of the SN ejecta (Dalcanton (2006)) 
(also called "the load factor" , adopted here as a free param- 
eter) and q is the yield. We fix the value of the yield to 
q — 0.02 and assume that the mass accreted to the disk has 

1 This is the escape velocity from R to infinity, not the escape 
velocity from the sites of star formation that are deeper in the po- 
tential well. This approximation does not affect the results for the 
massive spheroids since even the shallower potential at R is deep 
enough to prevent winds from being effective, see the discussion 
in section [5] For less massive systems we expect this approxima- 
tion to result in slightly smaller values for the parameter e than 
one would otherwise obtain. 



zero metallicity, i.e. we fix Zf = 0). Furthermore, we nor- 
malise the metallicities to the solar value for which we adopt 
Zq = 0.02. The chemical evolution of the ISM, described by 
Eq. (|18p contains three terms. The first one represents the 
chemical enrichment due to the evolution of stars. The sec- 
ond term represents the dilution of metallicity (if Zf < Z(t)) 
or the chemical enrichment (if Zf > Z(t)) of the ISM due 
to accreted material. The last term describes the dilution of 
metallicity (if r\ > 1) or the chemical enrichment (if 77 < 1) of 
the ISM due to galactic winds powered by stellar explosions. 

In recent theoretical work, what has been dubbed "the 
missing metals problem" has received considerable attention 
(see Prochaska et al. (2003) for an extensive discussion of 
this problem), namely the fact that the mean metallicity is 
~ 10 times lower than the value expected from the inferred 
star formation history. This problem may indicate a seri- 
ous flaw in our understanding of the interplay between star 
formation and metal enrichment. Therefore, we have intro- 
duced in our chemical evolution model an extra parameter 
/dii accounting phenomenologically for these effects. This al- 
lows the metallicity predictions of Eq. (|18p to be adjusted 
to match observational data. Thus we rescale the metallicity 
values given by solving Eq. (fTH|) by a factor f da , i.e. 

Z(t) = Z(t)/f da . (19) 

2.6 Summary of model parameters 

To summarize, our model is characterized by a set of 7 free 
parameters, that we denote by 8: 

9 = (log M min , e, x, v, t, 77, / di i) (20) 

The free parameters of the model (and the ones that we have 
chosen to fix) are summarized in Table [T] where we also give 
the prior ranges for our statistical analysis, i.e. the ranges 
within which their values are allowed to vary (see 14.21 for 
more details). 



3 INFLUENCE OF MODEL PARAMETERS ON 
OBSERVABLE QUANTITIES 

In this section, we discuss the impact of each of the 7 free 
parameters in our model (given in Eq. (|20p and in the up- 
per part of Table [1} on the physical observables introduced 
above, namely the SFR, SN type II rate, metallicity and 
baryonic fraction in structures. We also present a physical 
interpretation of the observed behaviour of these quantities. 
As a fiducial model we fix the parameter values to the follow- 
ing values: logM m i n = 8, e = 0.1, x = 1.7, r = 3, v = 1.4, 
rj — 10 and /dii = 2. We then proceed to vary one of the 
parameters at a time to get a feeling for the physical impact 
of each of them. 

Figure [1] shows the model dependence on the minimum 
mass of collapsed dark matter halos, logM m i n . Smaller val- 
ues of this parameter describe a scenario where star for- 
mation is hierarchical and follows the growth of structures, 
while higher values of log Mmtn correspond to star forma- 
tion occurring in massive spheroids. Correspondingly, for 
small logM m i n star formation begins earlier, as apparent 
from the top panel of Figure [l] At small redshift, a smaller 
log M m in leads to reduced SFR, since the relatively strong 
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Quantity 


Symbol 


Defined 


Prior range or value 


Minimum mass of collapsed haloes (-^^min i^ -^^"o) 


lop" M ■ 


Sec 


[2~2l 


5 < loe M ■ < 13 


SN type II energy efficiency factor 


e 


En 




0.01 < e ^ 0.45 


IMF power— law index 


X 


Rn 




3 < x ^ 2 


SFR normalization parameter (Gyr -1 ) 


V 


Rn 




o.oi s; v s; 5 


SFR timescale (Gyr) 


T 


En 




1 SC T ^ 5 


Winds load factor 


V 


Rn 




o s; t] s; 30 


Metals dilution factor 


./dil 


Rn 




1 /dil =S 30 


Baryon density parameter 


n b 






0.044 


Matter density parameter 








0.27 


Cosmological constant density parameter 


n A 






0.70 


Hubble constant (km/scc/Mpc) 








71 


Rms fluctuations amplitude 


cr.s 






0.9 


Dark matter to baryons bias parameter 


b 






1.0 


Minimum fraction of baryons when star formation begins 


/bar( 2 init) 






0.01 


Kinetic energy from stellar explosions 








10 51 ergs 


Yield 


<? 






0.02 


Metallicity of accreted material 


Z F 










Table 1. Upper part: free model parameters and priors used in the analysis. Top-hat (fiat) priors have been adopted on the parameter 
ranges indicated. Lower part: model parameters that have been fixed. 



winds (e = 0.1 in this example) drive the gas out of the sys- 
tem for shallower potentials. For large logM m i n , the build- 
up of metals is delayed in time but the metallicity can reach 
larger values, since supernova-powered winds are less impor- 
tant in massive systems (middle panel). As it is clear from 
Eq. ©, the percentage of baryons allocated to dark matter 
halos /bar increases for decreasing log M m i n (bottom panel). 
Due to the short lifetime of massive progenitors, the SN rate 
is essentially identical to the SFR, and we therefore do not 
display it. 

In Figure [5] we show the model sensitivity to the param- 
eter e, defining the percentage of supernovae energy that 
goes to the ISM. This parameter essentially describes the 
strength of galactic winds. The physical interpretation of 
high values of e is that strong winds driven by feedback en- 
ergy are maintained in dark matter halos. For increasing 
value of e, galactic winds become stronger and the star for- 
mation rate is reduced since less gas is available to make 
stars (top panel). This effect is more important for the shal- 
lower gravitational potential of the low mass halos, i.e. for 
smaller log A/ m i n (in this example, log Afrnin = 8) . As already 
remarked above, higher values of e corresponds to winds 
sweeping out metals from the ISM and thus to lower metal- 
licity (bottom panel). Again, this effect is most important 
for the low mass halos where their gravitational potential is 
relatively shallow. 

The sensitivity to the parameter x, giving the slope of 
the initial mass function, is shown in Figure [3] We can see a 
strong influence of x on the supernovae type II rates (middle 
panel), a consequence of Eq. (|17[l . Decreasing the value of 
x (i.e., making the IMF shallower) corresponds to a larger 
number of more massive stars, and hence the supernovae 
type II rate increases. Taking into account that in our model 
each supernova gives a constant percentage of its energy to 
the ISM, small values of x result in stronger galactic winds 
and so in smaller star formation rates (top panel) and a less 
enriched ISM, hence smaller metallicity (bottom panel). For 
the extreme case that x = 1 (very flat IMF) the supernovae 



type II rate is very large at high redshift causing very strong 
winds that reduce the SFR quickly. This causes the spike in 
Figure [3] (middle panel). 

Figure |3] shows the model sensitivity to the parameter 
p, connected with star formation through the proportional- 
ity factor that defines the efficiency of star formation, see 
Eq. (I16[l . Increasing the value of v the star formation be- 
comes more efficient and the interstellar medium becomes 
highly enriched in metals by evolving stars. On the contrary, 
smaller values of v lead to a less efficient star formation. 

The influence of the parameter r, defining the charac- 
teristic timescale of star formation, is displayed in Figure [5] 
Decreasing the value of r leads to the star formation activity 
ending sooner and to an ISM which is therefore poorer in 
metals. Larger values of r result in an enriched ISM since 
galaxies are active, in terms of star formation, for a longer 
period. 

The influence of the parameter 77, controlling the metal- 
licity of the ejecta, is displayed in Figure |6] A larger value 
of r\ leads to a decrease in metallicity of the system, since 
the metallicity of the winds is increased by a factor of r\ 
wrt the mean metallicity, see Eq. (JTHJ . Finally, we do not 
display the impact of the dilution factor /an, since its value 
merely rescales the metallicity by a multiplicative factor, see 
Eq. (JIl}. 

We now turn to discuss the data employed and the de- 
tails of our statistical treatment and fitting procedure. 



4 DATA AND STATISTICAL ANALYSIS 

Combining different types of observations to maximise their 
constraining power on multi-dimensional parameter spaces 
has become a common approach in cosmology. Following 
an approach similar in spirit, in this work we perform a 
simultaneous analysis of star formation history, SN rates, 
metallicity and baryonic fraction data in order to find tight 
constraints on the parameters of our model, Eq. (|20p . One 
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Figure 1. Dependence of the SFR, metallicity and baryonic frac- 
tion in structures (panels from top to bottom) on the minimum 
mass of collapsed dark matter halos, log M m j n . The curves are for 
log M min = 6, 8, 10, 12, from thin to thick. 



of the aim of this paper is to provide the first complete 
statistical analysis of existing metallicity, SFR, SNII rate 
and local collapse baryon fraction data in a realistic model. 
Wc first describe the data employed in section |4~T1 then we 
outline the Bayesian fitting procedure that vastly improves 
on usual fixed-grid scans in section |4~21 




0.1 1.0 10.0 

Redshift 

Figure 2. Dependence of the SFR (top panel) and of the 
metallicity (bottom panel) on the model parameter e, describ- 
ing the strength of galactic winds. The curves are for e = 
0.05,0.15,0.25,0.35, from thin to thick. 



4.1 Observational constraints 

The usual compilations of measurements for the SFR and 
metallicity observations (such as the ones used e.g. in Daigne 
et al. (2005) are unsuitable for a robust statistical analysis, 
because of the large systematic differences among measure- 
ments at about the same redshift performed over a range 
of different systems. In fact, when using such a "raw" data 
compilation the statistical fit is usually dominated by only 
a few data points with very small errorbars, while the large 
majority of observations carry almost no statistical weight. 
This is clearly less then satisfactory. To cure this effect, it be- 
comes important to bin the observations in such a way as to 
account for possible systematic uncertainties among differ- 
ent measurements at the same redshift. This problem is ad- 
dressed here for the first time by employing a Bayesian pro- 
cedure that accounts for possible systematic differences be- 
tween measurements, based on the treatment given in (Press 
(1996)). The details of the method are given in Appendix 1X1 
We apply the binning procedure described in the Ap- 
pendix to the data points for the metallicity in the inter- 
stellar medium (ISM) given by Prochaska et al. (2003). By 
using Eq. (|A3|I we place the 125 measurements in 8 bins, 
ranging in redshift from z — 0.85 to z = 4.45. The bin 
distribution and spacing has been chosen to obtain a rea- 
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Figure 3. Dependence of the SFR (top panel), SN rate (mid- 
dle panel) and of the metallicity (bottom panel) on the model 
parameter x, controlling the slope of the IMF. Curves are for 
x = 1.0, 1.4, 1.7, 2.0, from thin to thick. 



Figure 4. Dependence of the SFR (top panel) and the metallicity 
(bottom panel) on the parameter u, controlling the efficiency of 
star formation. The curves are for v = 0.5, 2, 3.5, 5, from thin to 
thick. 



Redshift 


Metallicity 


Number 




[M/H]/[M/H] Q 


of points 


0.85 


-0.83 ±0.11 


6 


1.45 


-1.06 ±0.09 


4 


1.95 


-0.93 ±0.12 


17 


2.45 


-1.36 ±0.14 


29 


2.95 


-1.56 ±0.19 


18 


3.45 


-1.78 ±0.08 


28 


3.95 


-1.80 ±0.06 


16 


4.45 


-1.76 ±0.11 


7 



Table 2. Binned measurements of [M/H] number density relative 
to solar values, after the statistical treatment of the data (see 
Appendix IA1I for details). 



sonable large number of points in each bin, while simulta- 
neously having a sufficiently small redshift spacing between 
bins. The measurements of [M/H] number density relative 
to solar metallicity obtained after the statistical rebinning, 
are summarized in Table [5] 

For the case of cosmic SFR data, our statistical rebin- 
ning is modified in order to take into account the redshift un- 
certainty in the raw data. Details are given in Appendix lA2l 
We take the compilation of "raw" data out to z ~ 5 from 



Hopkins (2004), excluding only one measurement corre- 
sponding to the cosmic star formation at z = 0.005T0.005, 
reported by Condon (1989). Instead, we replace this point 
by more recent measurement at the same redshift as re- 
ported by the same author (Condon et al. (2002)). Both 
these measurements use as cosmic star formation estimator 
counts at 1.4 GHz. From the raw data we derive binned val- 
ues in 12 redshift bins, with centers ranging fron z — 0.035 
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Figure 5. Sensitivity of the model to the parameter t, giving 
the characteristic timescale of star formation. The curves are for 
r = 1, 2, 3, 4, from thin to thick. 



to z = 5.12 by using Eq. l|A6f> . The resulting bins with their 
errors are summarized in Table [3] 

Furthermore, the SFR predictions of our model are cor- 
rected to account for dust absorption. There are large un- 
certainties associated with dust absorption correction, this 
is why at low redshift (i.e., for bins with Zb ^ 3) we employ 
both a "normal dust correction" of 1.0 mag and a "large 
dust correction" of 1.8 mag. These two choices are made in 
view of the fact that they seem to bracket the expected val- 
ues valid over a broad range of systems (Schiminovich at al. 
(2005)). For bins at a higher redshift (zt > 3) we adopt a 
fixed dust correction of 0.4 mag, following Schiminovich at 
al. (2005) . We shall see in the next section that the dust ab- 
sorption correction scheme one adopts has a crucial impact 
on the resulting physical scenario. 

The present-day fraction of baryons in structures, as 
estimated by Fukugita & Peebles (2004) is taken to be 



/bar(z = 0) = 0.61±0.11. 



(21) 



The data for the core collapse supernovae are 
taken from the Great Observatories Origins Deep Sur- 
vey (GOODS, Dahlen et al. (2004)). The GOODS core 
collapse supernovae rates have been placed in two bins at 
z = 0.3±0.2 and z = 0.7±0.2. For the local rate (at z = 0) 
we adopt the value from Cappellaro et al. (1999). We con- 
vert the local rate from supernovae units as described in 



Figure 6. Sensitivity of the metallicity to the parameter r\. 
Curves are for 77 = 5, 10, 15, 20, from thin to thick. 



Redshift 


SFR density 




logQaOIMsyr^Mpc- 3 ] 


0.012 


-1.78 ±0.18 


0.135 


-1.45 ±0.07 


0.275 


-1.45 ±0.06 


0.405 


-1.37 ±0.23 


0.580 


-1.08 ±0.08 


0.755 


-1.03 ±0.06 


0.905 


-0.98 ±0.07 


1.150 


-0.92 ±0.09 


1.650 


-0.63 ±0.26 


2.520 


-0.63 ±0.21 


3.770 


-0.79 ±0.10 


5.120 


-0.88 ±0.29 



Table 3. SFR density data after our statistical binning of the the 
"raw" SFR data compilation (see Appendix IA2I for details). No 
dust correction has been applied to this values. 



Dahlen et al. (2004). The 3 above mentioned data points 
are summarized in Table [4] 



4.2 Bayesian Markov Chain Monte Carlo analysis 

After the statistical rebinning of the data described above, 
the likelihood function P[d\ff) is the sum of four independent 
terms, describing the observations of the SFR, the metallic- 



Redshift SN type II rate 
z [Sne yr -1 Mpc -3 ] 



0.0 
0.3 
0.7 



6.16 ±2.92 
26. 20 4 



-,+7.83 
-9.18 



41.32; 



-11.06 
-10.75 



Table 4. Measurements of supernovae type II rate as a function 
of redshift. 
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ity, the SN rate and the baryonic fraction, respectively: 

P(d|(9) = £SFR +£mct +£sn +£(,- (22) 

We model each of the above four terms as a product of the 
data points for each observable, taken to be independent 
and with Gaussian noise 



the chain. This is achieved by defining the transition kernel 

as 



Xobe 



= -2 In Cobs = X> 



(23) 



where "obs" stand for SFR, metallicity, SN or baryon frac- 
tion in structures, and the means di and standard deviations 
Ui of the data points are given in Tables [2HH and Eq. (121 j) . 
The normalization constant does not matter, as we are only 
interested in the relative posterior probability density, as we 
now discuss. 

From the likelihood function of Eq. (|23p we obtain the 
posterior probability for the parameters of interest, P(6\d), 
via Bayes' theorem, 

P(d\6)P(6) 



P{0\d) = 



P(d) 



(24) 



where P{9) is the prior probability distribution ( "prior" for 
short) and P(d) is a normalization constant that does not 
depend on the parameters and can therefore be neglected 
in the following (see Trotta (2007a); Trotta (2007b) for 
more details on Bayesian parameter inference and model 
comparison). We adopt flat (i.e., top-hat) priors on our set 
of parameters 9 given in Eq. (|20p in the ranges given in Ta- 
ble [T] which means that the posterior probability distribu- 
tion function (pdf) is simply proportional to the likelihood. 

In order to explore efficiently our our 7-dimensional 
parameter space, we employ a Markov Chain Monte Carlo 
(MCMC) procedure, with some of the routines adapted from 
the publicly available COSMOMC packagfl The great advan- 
tages of MCMC methods are that the computational time 
scales approximately linearly with the number of dimensions 
of the parameter space, and that the marginalized posterior 
distribution for the parameters of interest and their correla- 
tions can be simply recovered by plotting histograms of the 
sample list. We follow the procedure outlined in de Austri 
et al. (2006), to which we refer for further details. Here we 
only briefly sketch the main points. 

The aim of an MCMC is to produce a series of sam- 
ples in parameter space (a Markov Chain) with the property 
that the density of points is proportional to the probability 
distribution (the target density) one is interested in map- 
ping, in our case the posterior pdf of Eq. (|24|l . There are 
several algorithms that can produce a chain with the re- 
quired properties. Here we employ the Metropolis-Hastings 
algorithm (Metropolis et al. (1953); Hastings (1970)): the 
chain is started from a random point in parameter space, 
6o, and a new point 9\ is proposed with an arbitrarily pro- 
posal density distribution q(9 n , (9 n +i). The transition kernel 
T(9 n ,9 n +i) gives the conditional probability for the chain 
to move from 9 n to # n +i, and it must satisfy the "detailed 
balance" condition 

P(6 n+ i\d)T(6 n+ i,9 n ) = P(6 n \d)T(e n ,e n +i) (25) 

so that the posterior P{0\d) is the stationary distribution of 

2 Available from cosmologist . info 



T(0 n ,0n+i) = q(9 n , 9 n +\)a(6n,0n+i), 



(26) 



a(9 n ,9 n+1 ) = mm jl, p^^^ ) , (27) 

where a(9 n , 9 n +i) gives the probability that the new point 
is accepted. Since P{0\d) oc P(d\9)P(9) and for the 
usual case of a symmetric proposal density, q(0 n ,9 n +i) = 
q(9 n +i,9 n ), the new step is always accepted if it improves 
on the posterior, otherwise it is accepted with probability 
P(d\9 n+ i)P{9 n+1 )/ P(d\9 n )P(9 n ). The result is a sample list 
from the target distribution, from which all the statistical 
quantities of interest can readily be evaluated. Further de- 
tails about MCMC methods can be found e.g. in MacKay 
(2003). 

Our Bayesian MCMC analysis allows us to not only to 
determine efficiently the best-fit value of the parameters, 
but also to explore correlations between the model parame- 
ters and estimate marginalized high probability regions, to 
which we now turn our attention. 



5 RESULTS AND DISCUSSION 

As mentioned above, we investigate two different dust cor- 
rection schemes for SFR data at low (z < 3) redshift, one 
termed "normal dust correction" and the other "high dust 
correction" . This is expected to roughly bracket the range of 
possible corrections. The outcome of our analysis is strongly 
dependent on which dust correction one chooses to employ, 
with the normal dust correction implying hierarchical star 
formation, while the high dust correction favours the mono- 
lithic scenario. 



5.1 Best fit models and parameter constraints 

The values of the best-fit model parameters for both dust 
correction schemes are given in Table[5] and the correspond- 
ing SFR, SN rate, metallicity evolution and baryonic frac- 
tion in structures are shown in Figure [7] The 1-dimensional 
posterior probability distributions (with all other parame- 
ters marginalized, i.e., integrated over) are plotted in Fig- 
ure [5] 

We first discuss the case with the normal dust correc- 
tion applied. In order to fit the (dust-corrected) SFR at both 
high and small redshifts, the model requires a small minimal 
mass (logM m i n ~ 6) and strong winds (e ~ 0.3). Although 
the value of the supernova energy transfer parameter is quite 
large, it is not too far away from theoretical predictions, 
which give an upper limit of e = 0.22 (Larson (1974)). An 
IMF power-law index x ~ 1.8, slightly larger then the Scalo 
IMF, is also preferred, which translates in fewer available 
supernovae. This is linked to the high value of e, since the 
energy transfer is so efficient that a large number of super- 
novae is not needed to get the appropriate feedback energy 
to reproduce the data sets. The metallicity load factor 77 can 
be connected with the IMF power-law index x and Dalcan- 
ton (2006) gives 77 values for a variety of IMFs. The value 
of r) for the Scalo (1986) IMF (x = 1.7, close to our best fit 
value, x — 1.77) is 77 = 16.8 — 18.6, in reasonable agreement 
with our value, rj = 8.74. This leaves metal-rich outflows as 
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Figure 7. Best-fitting models for the normal (solid line, hierarchical star formation, \ 2 = 26.60) and high (dashed, monolithic scenario, 
X 2 = 33.3) dust corrections, with parameters as in Table [3] In the top left panel, showing the SFR, the low redshift (z ^ 3) data have 
been corrected for dust employing a normal dust correction (1.0 mag, lower data points) or a high dust correction (1.8 mag, upper data 
points). The high dust correction data have bee shifted slightly to the right for display purposes. 



the only viable mechanism for producing the low effective 
yields observed in gas-rich galaxies, in agreement with sug- 
gestions presented in Dalcanton (2006). The dilution factor 
/dii is of order 2, which again is very reasonable, given the 
complex physics this parameter is supposed to summarize. 
The value of the chi-square of the best fit model in this case 
is 26.60 for 17 degrees of freedom, which suggests that our 
model captures the essential features of the data. Figure [7] 
shows the best fit models for normal dust correction (solid 
line) and high dust correction (dashed line). Both models 
provide an acceptable fit to the data, although in the nor- 
mal dust correction case the low-redshift metallicity and the 
present-day baryon fraction in structures appear in better 
agreement with the data. For a redshift above z ~ 5, the 
metallicity of the hierarchical model drops very sharply to 
because of the very significant winds. 

Turning now to the high dust correction case, we notice 
that the preferred values of the parameters in our model are 
very different from the previous case. Most importantly, a 
high dust correction at small redshift boosts the value of 
the SFR for z ^ 3, and this pushes our model to very large 
values of logA/ m i n , of the order logA/ m i n ~ 11 — 12. This 
implies that star formation occurs monolithically in heavy 



spheroids, as discussed in the introduction. We expect dry 
mergers to play a significant role in the build-up of mas- 
sive log M m in ~ 13 ellipticals in agreement with observa- 
tions showing that present-day spheroidal galaxies on aver- 
age have undergone between 0.5 and 2 major dry mergers 
since z ~ 0.7 (Bell et al. (2006)). Furthermore, we see from 
Figure[7] (dashed curves) that the onset of both the SFR and 
metals build-up is significantly delayed in this scenario, un- 
til about z ~ 5. The supernovae energy transfer parameter 
e becomes essentially irrelevant for such large values of the 
minimum mass, since the potential is deep enough to retain 
the ejected gas. The peak in the probability distribution for 
e observed in Figure [5] is therefore mostly a consequence of 
a volume effect of our Bayesian MCMC scanning technique. 
The star formation timescale r ~ 3.5 Gyr is in good agree- 
ment with theoretical models for Milky Way size disk galax- 
ies (with virial mass log M v i r close to our best fit value for 
logM m in). The IMF index is tilted towards extreme values, 
thus reducing the SN rates but boosting the SFR (compare 
Figure [3]). This in turns increases the metallicity, and a large 
dilution factor, fan ~ 20, is required to bring the predictions 
in line with observations. We notice that this agrees within 
a factor of two with the value already found in previous 
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Figure 8. 1— dimensional marginalised posterior probability distributions of the model parameters (normalized to their peak values). 
Solid histograms arc for the normal dust correction case (hierarchical scenario), dotted for the high dust correction (monolithic model). 



Parameter 


Normal dust correction 
Best fit 68% range 95% range 


Best fit 


High dust correction 

68% range 95% range 


log M min 


5.17 


< 6.11 


< 7.33 


11.60 


(11.27, 11.61) 


(10.97, 11.61) 


e 


0.32 


> 0.30 


> 0.21 


0.23 


< 0.17 


0.39 


X 


1.77 


(1.73,1.82) 


(1.68, 1.87) 


1.97 


> 1.95 


> 1.90 


v (Gyr- 1 ) 


4.15 


> 2.20 


> 1.01 


1.81 


(1.29, 1.90) 


(1.04,2.25) 


r (Gyr) 


3.59 


(3.48,4.22) 


(3.23,4.69) 


3.65 


(3.45,4.07) 


(3.24,4.46) 


V 


8.74 


(4.21, 15.70) 


(2.24,24.49) 


0.02 


< 3.80 


15.59 


/dil 


2.72 


(1.84,5.52) 


(1.27,9.34) 


20.75 


(15.65,22.34) 


(13.06,26.28) 


x 2 

X 2 /dof 


26.60 
1.6 








33.3 
2.0 





Table 5. Best— fit parameter values and marginalized 68% and 95% intervals for the normal (1.0 mag for z ^ 3) and high (1.8 mag for 
z ^ 3) dust corrections. For cases where only an upper or lower limit is found within our prior ranges, we give one— tail intervals. We also 
give the best-fit chi-square and the reduced chi-square, where the number of degrees of freedom (dof) is 17 (for 24 data points and 7 
free parameters). 



works on the metallicity of SN ejecta, which was of order 10. 
However, the extremely steep IMF that this model prefers 
(x ^ 2) appears to exclude the possibility that stellar ex- 
plosions are the main mechanism that drive galactic winds. 
This is reasonable, since supernova-driven gas flows cannot 
escape from massive galaxies' potential wells. A resolution to 
this wind dilemma could come from the hypothesis of super- 
massive black hole (SMBH) induced outflows (Silk (2005)). 
In fact, the very low value of the load factor (rj = 0.02) 
is consistent with this scenario since the SMBH undergoes 
most of its growth in the gas rich phase and its outflow ex- 
pels mostly unprocessed gas. Although our model does not 
include the physics of SMBHs, it is tempting to say that our 
best fit model suggests that SMBHs should play a key role 
in the evolution of massive spheroids. 

In general, we observe that the high dust correction 



case seems to stretch our model parameters to extreme val- 
ues, suggesting either a strong tension between datasets 
(mostly SFR and metallicity data) or a failure of the model 
to fully encapsulate all of the relevant physical processes. 
Even though with a reduced chi-square per dof of 2.0 this 
scenario is less favoured than the hierarchical star formation 
model discussed above, it appears that the monolothic for- 
mation model cannot be dismissed yet. It is interesting that 
our 7 parameter model is able to describe both cases, and 
that the SFR dust correction plays a major role in defining 
which scenario is preferred. 



5.2 Correlations among parameters 

We now turn to discuss the most relevant correlations among 
the model parameters in light of their physical interpretation 
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Figure 9. Contours enclosing joint 2D 68% and 95% regions, with all other parameters marginalized, for both the "high dust correction" 
(dashed, monolithic scenario) and the "normal dust correction" case (solid, hierarchical star formation). 



and of their impact on the observables, as shown in section[3] 
Figure [5] shows a selection of 2D joint posterior probability 
distributions for logMmin, e, x, v, r, and /an, thus giving 
complementary information to the ID distributions plotted 
in Figure [5] The contours enclose joint 2D 68% and 95% 
regions, with all other parameters marginalized, for both 
the "high dust correction" (dashed) and the "normal dust 
correction" case (solid). 

In the first panel on the left of Figure [9] showing the 
x — 6 plane, we observe a positive correlation between the 
IMF power-law index and the SN type II energy efficiency 
factor. This is expected, since an IMF with a higher power- 
law index produces less SNe, each of which has to con- 
tribute more energy, leading to higher values for the param- 
eter e (compare Figures [2] and panels showing the SFR 
and metallicity dependence). For the "high dust correction" 
case (dashed lines) this correlation is weaker, confirming our 
conclusion that SNe type II cannot drive the winds in the 
massive spheroids. The e — log M m i n plane shows that for 
structures of smaller mass ("normal dust correction" case, 
solid lines) the parameter e needs to be large, while for high 
mass structures (as preferred in the "high dust correction" 
case, dashed curves), e is essentially unconstrained, indicat- 
ing that SNe feedback is irrelevant for massive spheroids. 
The different physical processes taking place in small struc- 
tures and massive spheroids can be further investigated by 
looking at the correlations in the fdu — v plane. We expect 
to find a positive correlation among v and the dilution fac- 
tor /dii, as larger v increases the SFR (compare Figure [4} 
thus leading to a more metal-rich ISM. To bring this back 
in line with the data, a larger dilution factor is needed. The 
above line of reasoning explains the strong positive corre- 
lation one observes for the high mass structures (dashed) 
where winds do not play a strong role and metals cannot 
escape from the structure. In contrast, metal-rich winds are 
dominant for smaller structures (solid curves), thus expelling 
most of the metals produced. This results in almost no cor- 
relation between /an and u, since the impact of u on the 
SFR and metallicity predictions can be mimicked by a dif- 
ferent combination of values for e and logM m i n . Finally, in 
the right-most panel of Figure [9] we display the probability 
distribution in the r — v plane, which exhibits a strong nega- 
tive correlation. Again, this is expected on the grounds that 
large values of the parameter v increase the SFR (compare 
Figure [4| and a smaller time-scale is thus required in order 
to quench star formation fast enough (see Figure [SJ). 



6 CONCLUSIONS 

We have presented a well-motivated physical model of the 
cosmic star formation incorporating supernova feedback, gas 
accretion and enriched outflows. We computed the cosmic 
star formation history and the chemical evolution in the in- 
terstellar medium of forming galaxies as a function of red- 
shift, and we presented for the first time a full statistical 
treatment of the observational data, which accounts for the 
possibility of systematic errors in the data sets. 

We have employed four different observational data — 
the observed cosmic star formation rate up to z ~ 5, the ob- 
served rate of type II supernova up to z ~ 0.7, the present 
baryon fraction in structures and the evolution of the metal 
content in the ISM — to derive constraints on the free pa- 
rameters of our model. After employing a Bayesian proce- 
dure to rebin the SFR and SN rate data, we found that 
the low redshift (z ^ 3) SFR dust correction adopted has a 
critical impact on the scenario favoured by the data. 

For what we have termed "normal dust correction", 
the hierarchical star formation model is preferred, where 
star formation occurs in small structures first and super- 
novae winds are important. While the winds load factor re- 
mains poorly constrained, we can conclude that larger values 
are preferred, in agreement with previous work (Dalcanton 
(2006)). Applying a larger dust correction at small redshifts, 
we found that the data on the contrary favour high val- 
ues for the minimum mass of a dark halo of the collapsed 
structures (monolithic star formation scenario). This case 
requires a large dilution factor, a rather extreme IMF slope 
and a fairly small winds load factor, at the model param- 
eters are pushed at the boundaries of the available range. 
We have suggested that this might be interpreted in terms 
of the presence of outflow from supermassive black holes, 
but this possibility will require further investigation. It is 
worth noticing that the monolithic star formation scenario 
has very little star formation beyond z ~ 5. Observations of 
the E-mode polarization power spectrum of the cosmic mi- 
crowave background however indicate that the Universe was 
re-ionized around z ~ 11 (Spergel et al. (2006)). This means 
that in this scenario the reionization mechanism has to be 
found elsewhere than in massive UV — emitting stars. Sev- 
eral alternatives have been explored in the literature, for ex- 
ample reionization by decaying particles (Hansen & Haiman 
(2004)), or a high-redshift population of mini-quasars that 
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can reionize the IGM up to 50% ionisation fraction (Dijkstra 
et al. (2004)). 

For both models, the IMF slope is large. Unfortunately, 
this does not help in distinguishing one model from the 
other, since observations have so far not yielded convincing 
results concerning the form of the stellar IMF or its varia- 
tions in space and time (Scalo (1998)). The most important 
difference among the IMFs is that the fraction of high-mass 
stars is larger for a shallower IMF. Since only high-mass 
stars emit significant amount of ultraviolet light, this results 
in a spectrum which is more shifted towards the ultraviolet 
for a typical galaxy with an e.g. Salpeter IMF (x — 1.35) as 
compared with a Scalo IMF (x — 1.7). In turn, this leads to a 
different reionization history, which can be in principle com- 
pared with the optical depth to reionization as inferred from 
cosmic microwave background polarization measurements. 

While the monolithic scenario is less preferred in terms 
of quality of fit, it is clear that more work is required to 
be able to draw firm conclusions as to the viability of the 
two different models. Of particular importance remains the 
statistical treatment of the data, for which we have here 
presented a new procedure that we hope will prove useful 
for future work. 
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APPENDIX A: BINNING OF DATA 
ACCOUNTING FOR UNDETECTED 
SYSTEM ATICS 

Al No redshift uncertainty 

We wish to define B bins in redshift space. Within each bin 
b, 1 ^ b ^ B, we have a collection of measurements (in our 
case, metallicity or SFR data), each with its own statisti- 
cal accuracy and possibly an unspecified systematic error. 
That systematic differences above the quoted statistical er- 
rors dominate the raw data is apparent from a plot of the 
unbinned metallicity or SFR observations, that show a scat- 
ter of up to an order of magnitude for observations at about 
the same redshift. The origin of the systematic discrepancy 
can vary, from underestimated statistical errors in the ob- 
servation to intrinsic dispersion in the observed systems to 
differences in the way the data are collected. In the presence 
of systematic errors, we cannot simply take the weighted 
average of the data within each bin. Instead, we model the 
presence of unknown systematics as follows. 
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Let us consider the measurement of a quantity yt in a 
top-hat bin b, 1 ^ b ^ B — in our case, this represents 
the metallicity value at the redshift of the bin, Zb, and we 
assume we can neglect the redshift uncertainty of the mea- 
surements (this issue is addressed in the next section). Each 
measurement consists of a central value di and a statistical 
error a, 1 ^ i ^ JVj, for iVj, different measurements within 
bin b. If the i-th datum does not suffer from a systematic er- 
ror (or where the systematic error, Si, is negligible compared 
with the quoted statistical error), the likelihood function is 
modelled as a Gaussian with the quoted standard deviation 



Pi,g{di\yb) = 



27TO". 



exp 



1 ^ di — y b 
2 



(Al) 



For the sake of brevity, let us denote such measurements 
as "good" measurements, as indicated by the subscript g. 
If the datum suffers from an undetected systematic, i.e. the 
dominant error is Si S> en, the likelihood is instead given by 
(neglecting the statistical error wrt the systematic one) 



Pi,s{di\yb) 



2nSi 



exp 



di — yb 
Si 



(A2) 



where the subscript s denotes "systematics" , or "spoiled" 
measurements, for brevity. Now of course we do not know 
which measurements suffer from systematic, but this can 
be determined statistically using the following procedure 
(adapted from Press (1996)). 

We denote by p the probability that each of the mea- 
surements i in bin b is a "good" one. Conversely, 1 — p is the 
probability that the datum suffers from systematics. Fur- 
thermore, we include a binary vector V = (Vi, . . . , Vjv b ), 
whose elements Vi (1 ^ i ^ JV b ) can either be or 1, deter- 
mining whether the datum i is a good one (for Vi — 1) or a 
spoiled one (for Vi — 0). We can then compute the posterior 
probability for the value of the observed quantity yb in bin b 
by multiplying the individual contributions of the measure- 
ments in the same bin and marginalizing over the unknown 
quantities p and V (see Eq. (16) in Press (1996)): 



P(y b \d b ) oc 



j d P n\pp g ,i+(i-p)p s , 



(A3) 



where d b denotes the collection of measurements in bin d, i.e. 
db = (di, • • • , djv b )- For the prior probability on p we have 
assumed a flat prior distribution between ^ p ^ 1 and 
the proportionality factor might be determined by requiring 
that the likelihood be normalized to unity, but this is not 
necessary in our application. The precise numerical value of 
the error associated with systematics, Si does not matter, as 
long as Si S> Oi. In our case, we take S% to be unity on a log 
scale, corresponding to one order of magnitude uncertainty 
on the observable. 

From the posterior distribution (|A3|I . the central value 
of the bin b is obtained as the peak of the distribution, 
while the standard deviation is defined as the range enclos- 
ing 68.4% (lcr range) of the probability. These values are 
given in Table [2] for the metallicity data, and are then used 
for the likelihood function employed in the fit of the model. 
Of course one could as well employ the full probability dis- 
tribution of Eq. (|A3[l as the likelihood function, but for sim- 
plicity we have summarized it as a Gaussian with mean and 
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Figure A2. Raw SFR data and the binned values after the sta- 
tistical treatment including redshift uncertainties. No dust cor- 
rection has been applied to the data at this stage. 



standard deviation computed as described above. The col- 
lection of raw, unbinned metallicity data and the resulting 
bins form our statistical treatment are shown in Figure lATI 



A2 Accounting for redshift uncertainty 

When the observed quantity suffers from a substantial red- 
shift uncertainty, as in the case of the SFR data, we need 
to take into account the redshift error in our binning pro- 
cedure, as this introduces a further uncertainty as to which 
bin a given datum belongs to. The above procedure is then 
modified as follows. 

The probability that an observations with central red- 
shift Zi and redshift uncertainty n belongs to the fe-th red- 
shift bin (centered at redshift Zb) is modelled as a Gaussian, 
i.e. 



P(zi\z b ) = 



2%Ti 



exp 



z b 



(A4) 



Given the uncertainty on the location of the measurements 
in redshift, it is now impossible to assign data points to 
top-hat bins. Instead, one needs to marginalize over all pos- 
sible assignments of data points among redshift bins, with 
each point's contribution weighted by the conditional prob- 
ability of Eq. (|A4[1 . For each bin, let us introduce a new 
binary variable, Z = (Zi, . . . , Zti), whose elements indicate 
whether the z-th datum (1 ^ i ^ N) belongs to the bin 
under consideration (Zi — 1) or not (Zi = 0). If we knew 
which datum belongs to which redshift bin, then we could 
assign an exact binary sequence to Z (this corresponds to 
the case considered in the previous section). Instead, we sum 
(marginalize) over all possibilities, writing for the posterior 
probability of the SFR value yb at redshift Zb, given d, the 
full collection of data points at all redshifts 

P{Vb\d) = J2 P ^»' Z \ d )=J2 P (Vb\ Z , d)P(z l \z b , Zi = 1), 
z z 

(A5) 

where the conditional probability P(yb\Z,d) is given by 
Eq. (IA3|) . given a specific assignment for Z. The sum over 
Z can be replaced by a product of binomial terms, so that 



© 0000 RAS, MNRAS 000, 000-000 



we finally obtain, using Eqs. (EH} and flUJ 

P(yt\d)<x / dp ]J{[pP s>i + (l-p)P a , i ]P(a i |z 6 )}-l. 

^ i=l 

(A6) 

Notice that the product is here over all points in the dataset, 
not just over the ones in a bin, as in Eq. (|A3[) . 

Since we include the full dataset for each bin, the result- 
ing errors are in principle correlated across bins. However, 
the Gaussian term of Eq. (|A4[1 ensures that only "nearby" 
points give a non-negligible contribution to the value of bin 
b. We therefore consider it acceptable to ignore the correla- 
tion among bins when using the mean and standard devia- 
tion of Eq. (1A6[) for the likelihood function for the SFR. The 
results from this procedure are tabulated in Table [3] and are 
plotted alongside with the raw, unbinned data in Figure lA"2l 



Monolithic or hierarchical star formation 
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